#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <memory.h>

//#define DBG_PRINT(_X_)  printf _X_
#define DBG_PRINT(_X_)

void gaus(double a[4][4],double *b, double *res)
{
	double s,c[4],x[4],wer;
	
	int i,j,l;
	DBG_PRINT(("\nOutput isxodnik\n"));

	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			DBG_PRINT(("\t%2.2f ",a[i][j]));
		};
		DBG_PRINT(("\t| %2.2f",b[i]));
		DBG_PRINT(("\n"));
	};

	for(l=0;l<4;l++)
	{
		for(i=l;i<4;i++)
		{
			c[i]=a[i+1][l]/a[l][l];
			b[i+1]=b[i+1]-b[l]*c[i];
			for(j=l;j<=3;j++)
			{
				a[i+1][j]=a[i+1][j]-a[l][j]*c[i];
			};
		};
	};

	DBG_PRINT(("\nTreugolnik:\n"));

	for(i=0;i<4;i++)
    {
		for(j=0;j<4;j++)
		{
			DBG_PRINT((" %2.2f ",a[i][j]));
		};
		DBG_PRINT((" | %2.2f",b[i]));
		DBG_PRINT(("\n"));
    };

	DBG_PRINT(("\n\n"));

	for(i=0;i<4;i++)     
	{
		x[i]=0;
	};

	for(i=3;i>=0;i--)
	{
 		wer=0;
		if(i!=4) for(j=i+1;j<4;j++)
 		{
 			wer=wer+a[i][j]*x[j];
 		};

 		x[i]=(b[i]-wer)/a[i][i];
 	};

 	DBG_PRINT(("Reshenie:\n"));
	for(i=0;i<4;i++)
		{
		DBG_PRINT(("x[%d] = %2.2f\n",i,x[i]));
		res[i]=x[i];
		};
}

void main()
{
	double b[]={0.11,2.00,0.11,0.13};
	double f[]={0.11,2.00,0.11,0.13};
	int i,j;
	double epsilon = 0.000001;
	double a[4][4]=
		{
			{0.17,0.75,-0.18,0.21},
			{0.75,0.13,0.11,1.00},
			{-0.33,0.11,3.01,-2.01},
			{0.11,1.12,1.11,-1.31}
		};
	double temp_f[4],f0[4],r0[4],xx[4],dx[4],xx0[4],delta_xx[4];
	bool isExit;

	for(i=0;i<4;i++)
    {
		for(j=0;j<4;j++)
		{
			printf(" %2.2f ",a[i][j]);
		};
		printf(" | %2.2f",b[i]);
		printf("\n");
    };

	for(i=0;i<4;i++)
		xx0[i]=0;

	//----------------------
	for(i=0;i<4;i++)
		temp_f[i]=f[i];
	gaus(a,temp_f,xx);


	int count = 0;
	do
	{
		DBG_PRINT(("\nf0:\n"));
		for(i=0;i<4;i++)		
		{		
			f0[i]=0;
			for(j=0;j<4;j++)
			{
				f0[i]=f0[i]+(a[i][j]*xx[j]); 
			};
			DBG_PRINT(("%2.2f ",f0[i]));
		};

		DBG_PRINT(("\nb:\n"));

		for(i=0;i<4;i++) {
			b[i]=(f0[i]-f[i])*(-1);
			DBG_PRINT(("%2.2f ",b[i]));
		}

		gaus(a,b,delta_xx);

		isExit = true;
		DBG_PRINT(("\ndelta:\n"));
		for(i=0;i<4;i++){
			xx0[i]=xx[i];
			xx[i]=xx[i]+delta_xx[i];
			DBG_PRINT(("%2.2f ",delta_xx[i]));
			if (delta_xx[i]>epsilon)
				isExit = false;
		}
		
		if (isExit == true)
			break;

		for(i=0;i<4;i++)
			f[i]=f0[i];
		count++;
	}
	while(count<200000);

	printf("\nPROVERKA______!!!!!!!!!!!!: count[%d]\n",count);
	for(i=0;i<4;i++)		
	{			
		f0[i]=0;
		for(j=0;j<4;j++)
		{
			f0[i]=f0[i]+(a[i][j]*xx[j]); 
		};
		printf("%2.2f ",f0[i]);
	};

	getch();
}